投票率55.93%はいかがですか？

Primitive Analysis

作者

所属

Jaehyun Song

Kansai University

1 セットアップ

コード

pacman::p_load(tidyverse, readxl, naniar, marginaleffects, modelsummary)

raw_df <- read_xlsx("Data/Data.xlsx")

df <- raw_df %>%
  replace_with_na(replace = list(Q15.1 = c(88, 99),
                                 Q15.2 = c(88, 99),
                                 Q16.1 = c(88, 99),
                                 Q16.2 = c(88, 99),
                                 Q2.1  = 99,
                                 Q2.3  = 99,
                                 Q2.4  = 99,
                                 Q3_1  = 99,
                                 Q4.1  = 99,
                                 Q7.1  = 99,
                                 Q13.1 = c(6, 99))) %>%
  mutate(Q15.1 = 6 - Q15.1,
         Q15.2 = recode(Q15.2,
                        "1" = "5",
                        "2" = "4",
                        "3" = "2",
                        "4" = "1",
                        "5" = "3"),
         Q16.1 = 6 - Q16.1,
         Q16.2 = recode(Q16.2,
                        "1" = "5",
                        "2" = "4",
                        "3" = "2",
                        "4" = "1",
                        "5" = "3"),
         Female = sex - 1,
         Age    = age,
         Educ   = Q2.1,
         Income = Q2.3,
         Urban  = 5 - Q2.4,
         Duty   = Q3_1,
         Interest = 5 - Q4.1,
         Independent = if_else(Q7.1 == 11, 1, 0),
         Voted       = if_else(Q13.1 <= 3, 1, 0)) %>%
  rename(Treatment = votingrate) %>%
  filter(Q3_8 == 4)

df2 <- df %>%
  group_by(Treatment) %>%
  summarise(Outcome = mean(Q16.1, na.rm = TRUE),
            .groups = "drop")

2 55.93%に対する評価

コード

fit1 <- lm(Q15.1 ~ Female + Age + Educ + Income + Urban + 
             Duty + Interest + Independent + Voted,
   data = df)

modelsummary(fit1,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.986 (0.174) ***
Female	0.141 (0.049) **
Age	-0.003 (0.002)
Educ	0.000 (0.026)
Income	0.004 (0.005)
Urban	0.001 (0.026)
Duty	-0.146 (0.022) ***
Interest	-0.133 (0.038) ***
Independent	-0.089 (0.049) +
Voted	-0.247 (0.071) ***
Num.Obs.	1723
R2	0.093
R2 Adj.	0.088
AIC	4746.3
BIC	4806.2
Log.Lik.	-2362.136
RMSE	0.95

コード

tidy(fit1, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(term = fct_inorder(term),
         term = fct_rev(term),
         sig  = if_else(p.value < 0.05, "sig", "insig")) %>%
  ggplot() +
  geom_vline(xintercept = 0) +
  geom_pointrange(aes(x = estimate, y = term,
                      xmin = conf.low, xmax = conf.high, color = sig)) +
  scale_color_manual(values = c("sig" = "black", "insig" = "gray70")) +
  labs(x = "係数（正の場合、ポジティブ評価）") +
  guides(color = "none") +
  theme_bw(base_size = 12) +
  theme(axis.title.y = element_blank())

3 実験（処置効果）

コード

fit2 <- lm(Q16.1 ~ poly(Treatment, 2) + 
             Female + Age + Educ + Income + Urban + 
             Duty + Interest + Independent + Voted, 
   data = df)

modelsummary(fit2,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.399 (0.159) ***
poly(Treatment, 2)1	42.939 (1.023) ***
poly(Treatment, 2)2	16.364 (1.024) ***
Female	0.111 (0.045) *
Age	0.001 (0.001)
Educ	0.014 (0.025)
Income	-0.003 (0.004)
Urban	0.006 (0.024)
Duty	-0.104 (0.021) ***
Interest	-0.034 (0.035)
Independent	-0.021 (0.045)
Voted	-0.015 (0.065)
Num.Obs.	1761
R2	0.544
R2 Adj.	0.541
AIC	4629.5
BIC	4700.6
Log.Lik.	-2301.736
RMSE	0.89

コード

predictions(fit2, newdata = datagrid(Treatment = 20:80)) %>% 
  ggplot(aes(x = Treatment, y = estimate)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_point(data = df2, aes(x = Treatment, y = Outcome),
             color = "gray80", size = 3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI",
       caption = "注：点は投票率ごとの応答変数の平均値を示す。") +
  theme_bw(base_size = 12)

4 条件付き処置効果

4.1 性別

コード

fit3 <- lm(Q16.1 ~ poly(Treatment, 2) * Female + 
             Age + Educ + Income + Urban + 
             Duty + Interest + Independent + Voted, 
   data = df)

modelsummary(fit3,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.396 (0.159) ***
poly(Treatment, 2)1	41.346 (1.439) ***
poly(Treatment, 2)2	17.152 (1.441) ***
Female	0.112 (0.045) *
Age	0.001 (0.001)
Educ	0.014 (0.025)
Income	-0.003 (0.004)
Urban	0.007 (0.024)
Duty	-0.104 (0.021) ***
Interest	-0.032 (0.035)
Independent	-0.022 (0.045)
Voted	-0.018 (0.065)
poly(Treatment, 2)1 × Female	3.217 (2.049)
poly(Treatment, 2)2 × Female	-1.580 (2.052)
Num.Obs.	1761
R2	0.545
R2 Adj.	0.541
AIC	4630.4
BIC	4712.5
Log.Lik.	-2300.204
RMSE	0.89

予測値

コード

predictions(fit3,
            newdata = datagrid(Treatment = 20:80,
                               Female     = 0:1)) %>%
  mutate(Female = factor(Female, levels = 0:1, 
                         labels = c("男性", "女性"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Female), alpha = 0.25) +
  geom_line(aes(color = Female), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "性別", fill = "性別") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit3, variables = "Treatment",
                newdata = datagrid(Female = 0:1)) %>%
  mutate(Female = factor(Female, levels = 0:1, 
                         labels = c("男性", "女性"))) %>%
  ggplot(aes(x = Female, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "性別", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit3, variables = "Female",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "性別の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.2 年齢

コード

fit4 <- lm(Q16.1 ~ poly(Treatment, 2) * Age + 
             Female + Educ + Income + Urban + 
             Duty + Interest + Independent + Voted, 
   data = df)

modelsummary(fit4,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.416 (0.158) ***
poly(Treatment, 2)1	31.986 (3.388) ***
poly(Treatment, 2)2	6.091 (3.351) +
Age	0.001 (0.001)
Female	0.110 (0.045) *
Educ	0.014 (0.024)
Income	-0.003 (0.004)
Urban	0.006 (0.024)
Duty	-0.107 (0.020) ***
Interest	-0.035 (0.035)
Independent	-0.022 (0.045)
Voted	-0.024 (0.065)
poly(Treatment, 2)1 × Age	0.216 (0.063) ***
poly(Treatment, 2)2 × Age	0.204 (0.063) **
Num.Obs.	1761
R2	0.550
R2 Adj.	0.546
AIC	4611.3
BIC	4693.5
Log.Lik.	-2290.673
RMSE	0.89

予測値

コード

predictions(fit4,
            newdata = datagrid(Treatment = 20:80,
                               Age       = seq(20, 80, by = 20))) %>%
  mutate(Age = factor(Age, levels = seq(20, 80, by = 20), 
                      labels = paste0(seq(20, 80, by = 20), "歳"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Age), alpha = 0.25) +
  geom_line(aes(color = Age), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "年齢", fill = "年齢") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit4, variables = "Treatment",
                newdata = datagrid(Age = 18:79)) %>%
  ggplot(aes(x = Age, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "年齢（歳）", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit4, variables = "Age",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "年齢の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.3 教育水準

コード

fit5 <- lm(Q16.1 ~ poly(Treatment, 2) * Educ + 
             Female + Age + Income + Urban + 
             Duty + Interest + Independent + Voted, 
   data = df)

modelsummary(fit5,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.401 (0.159) ***
poly(Treatment, 2)1	42.410 (3.685) ***
poly(Treatment, 2)2	19.000 (3.745) ***
Educ	0.014 (0.025)
Female	0.110 (0.045) *
Age	0.001 (0.001)
Income	-0.003 (0.004)
Urban	0.006 (0.024)
Duty	-0.104 (0.021) ***
Interest	-0.035 (0.035)
Independent	-0.021 (0.045)
Voted	-0.014 (0.065)
poly(Treatment, 2)1 × Educ	0.159 (1.069)
poly(Treatment, 2)2 × Educ	-0.785 (1.073)
Num.Obs.	1761
R2	0.544
R2 Adj.	0.541
AIC	4632.9
BIC	4715.0
Log.Lik.	-2301.453
RMSE	0.89

予測値

コード

predictions(fit5,
            newdata = datagrid(Treatment = 20:80,
                               Educ      = 1:5)) %>%
  mutate(Educ = factor(Educ, levels = 1:5, 
                       labels = c("中学校", "高校",
                                  "短大・高専・専門学校",
                                  "大学", "大学院"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Educ), alpha = 0.25) +
  geom_line(aes(color = Educ), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "最終学歴（在学中を含む）", 
       fill = "最終学歴（在学中を含む）") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit5, variables = "Treatment",
                newdata = datagrid(Educ = 1:5)) %>%
  ggplot(aes(x = Educ, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  scale_x_continuous(breaks = 1:5,
                     labels = c("中学校", "高校",
                                "短大・高専・専門学校",
                                "大学", "大学院")) +
  labs(x = "最終学歴（在学中を含む）", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit5, variables = "Educ",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "最終学歴の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.4 世帯収入

コード

fit6 <- lm(Q16.1 ~ poly(Treatment, 2) * Income + 
             Female + Age + Educ + Urban + 
             Duty + Interest + Independent + Voted, 
   data = df)

modelsummary(fit6,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.399 (0.159) ***
poly(Treatment, 2)1	43.191 (1.615) ***
poly(Treatment, 2)2	17.810 (1.615) ***
Income	-0.003 (0.004)
Female	0.112 (0.045) *
Age	0.001 (0.001)
Educ	0.015 (0.025)
Urban	0.006 (0.024)
Duty	-0.104 (0.021) ***
Interest	-0.034 (0.035)
Independent	-0.022 (0.045)
Voted	-0.012 (0.065)
poly(Treatment, 2)1 × Income	-0.034 (0.193)
poly(Treatment, 2)2 × Income	-0.222 (0.193)
Num.Obs.	1761
R2	0.544
R2 Adj.	0.541
AIC	4632.1
BIC	4714.2
Log.Lik.	-2301.037
RMSE	0.89

予測値

コード

predictions(fit6,
            newdata = datagrid(Treatment = 20:80,
                               Income    = c(1, 10, 21))) %>%
  mutate(Income = factor(Income, levels = c(1, 10, 21))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Income), alpha = 0.25) +
  geom_line(aes(color = Income), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "世帯収入（万円）", fill = "世帯収入（万円）") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit6, variables = "Treatment",
                newdata = datagrid(Income = 1:21)) %>%
  ggplot(aes(x = Income, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "世帯収入", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit6, variables = "Income",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "世帯収入の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.5 居住地域

コード

fit7 <- lm(Q16.1 ~ poly(Treatment, 2) * Urban + 
             Female + Age + Educ + Income + 
             Duty + Interest + Independent + Voted, 
   data = df)

modelsummary(fit7,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.398 (0.159) ***
poly(Treatment, 2)1	38.079 (3.507) ***
poly(Treatment, 2)2	14.729 (3.467) ***
Urban	0.007 (0.024)
Female	0.112 (0.045) *
Age	0.001 (0.001)
Educ	0.014 (0.025)
Income	-0.003 (0.004)
Duty	-0.104 (0.021) ***
Interest	-0.035 (0.035)
Independent	-0.023 (0.045)
Voted	-0.015 (0.065)
poly(Treatment, 2)1 × Urban	1.622 (1.122)
poly(Treatment, 2)2 × Urban	0.539 (1.114)
Num.Obs.	1761
R2	0.545
R2 Adj.	0.541
AIC	4631.1
BIC	4713.3
Log.Lik.	-2300.573
RMSE	0.89

予測値

コード

predictions(fit7,
            newdata = datagrid(Treatment = 20:80,
                               Urban     = 1:4)) %>%
  mutate(Urban = factor(Urban, levels = 1:4, 
                        labels = c("町村", "人口10万人未満の市",
                                   "人口10万人以上の市",
                                   "政令市・23区"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Urban), alpha = 0.25) +
  geom_line(aes(color = Urban), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "居住地域", fill = "居住地域") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit7, variables = "Treatment",
                newdata = datagrid(Urban = 1:4)) %>%
  ggplot(aes(x = Urban, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  scale_x_continuous(breaks = 1:4,
                     labels = c("町村", "人口10万人未満の市",
                                "人口10万人以上の市", "政令市・23区")) +
  labs(x = "最終学歴（在学中を含む）", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit7, variables = "Urban",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "居住地域の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.6 投票義務感

コード

fit8 <- lm(Q16.1 ~ poly(Treatment, 2) * Duty + 
             Female + Age + Educ + Income + Urban + 
             Interest + Independent + Voted, 
   data = df)

modelsummary(fit8,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.402 (0.159) ***
poly(Treatment, 2)1	45.485 (3.337) ***
poly(Treatment, 2)2	12.301 (3.313) ***
Duty	-0.104 (0.021) ***
Female	0.112 (0.045) *
Age	0.001 (0.001)
Educ	0.014 (0.025)
Income	-0.003 (0.004)
Urban	0.006 (0.024)
Interest	-0.033 (0.035)
Independent	-0.022 (0.045)
Voted	-0.018 (0.065)
poly(Treatment, 2)1 × Duty	-0.672 (0.830)
poly(Treatment, 2)2 × Duty	1.070 (0.827)
Num.Obs.	1761
R2	0.545
R2 Adj.	0.541
AIC	4631.3
BIC	4713.4
Log.Lik.	-2300.626
RMSE	0.89

予測値

コード

predictions(fit8,
            newdata = datagrid(Treatment = 20:80,
                               Duty      = c(1, 3, 5))) %>%
  mutate(Duty = factor(Duty, levels = c(1, 3, 5), 
                       labels = c("そう思う", "どちらともいえない",
                                  "まったくそう思わない"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Duty), alpha = 0.25) +
  geom_line(aes(color = Duty), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "自分一人くらい投票しなくてもかまわない", 
       fill = "自分一人くらい投票しなくてもかまわない") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit8, variables = "Treatment",
                newdata = datagrid(Duty = 1:5)) %>%
  mutate(Duty = factor(Duty, levels = 1:5, 
                       labels = c("そう思う", 
                                  "ややそう思う",
                                  "どちらともいえない",
                                  "あまりそう思わない",
                                  "まったくそう思わない"))) %>%
  ggplot(aes(x = Duty, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "選挙では大勢が投票するものだから、自分一人くらい投票しなくてもかまわない", 
       y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit8, variables = "Duty",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "投票義務感の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.7 政治関心

コード

fit9 <- lm(Q16.1 ~ poly(Treatment, 2) * Interest + 
             Female + Age + Educ + Income + Urban + 
             Duty + Independent + Voted, 
   data = df)

modelsummary(fit9,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.405 (0.159) ***
poly(Treatment, 2)1	44.114 (4.195) ***
poly(Treatment, 2)2	6.940 (4.335)
Interest	-0.031 (0.035)
Female	0.110 (0.045) *
Age	0.001 (0.001)
Educ	0.015 (0.025)
Income	-0.003 (0.004)
Urban	0.005 (0.024)
Duty	-0.104 (0.021) ***
Independent	-0.018 (0.045)
Voted	-0.021 (0.065)
poly(Treatment, 2)1 × Interest	-0.391 (1.448)
poly(Treatment, 2)2 × Interest	3.331 (1.490) *
Num.Obs.	1761
R2	0.545
R2 Adj.	0.542
AIC	4628.4
BIC	4710.5
Log.Lik.	-2299.194
RMSE	0.89

予測値

コード

predictions(fit9,
            newdata = datagrid(Treatment = 20:80,
                               Interest  = c(1, 4))) %>%
  mutate(Interest = factor(Interest, levels = c(1, 4), 
                           labels = c("まったく関心がない", 
                                      "とても関心がある"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Interest), alpha = 0.25) +
  geom_line(aes(color = Interest), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "政治関心", fill = "政治関心") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit9, variables = "Treatment",
                newdata = datagrid(Interest = 1:4)) %>%
  mutate(Interest = factor(Interest, levels = 1:4, 
                           labels = c("まったく関心がない", 
                                      "あまり関心がない",
                                      "ある程度関心がある",
                                      "とても関心がある"))) %>%
  ggplot(aes(x = Interest, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "政治関心", 
       y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit9, variables = "Interest",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "政治関心の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.8 支持政党

コード

fit10 <- lm(Q16.1 ~ poly(Treatment, 2) * Independent + 
             Female + Age + Educ + Income + Urban + 
             Duty + Interest + Voted, 
   data = df)

modelsummary(fit10,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.396 (0.159) ***
poly(Treatment, 2)1	41.109 (1.334) ***
poly(Treatment, 2)2	16.945 (1.326) ***
Independent	-0.020 (0.045)
Female	0.109 (0.045) *
Age	0.001 (0.001)
Educ	0.015 (0.025)
Income	-0.003 (0.004)
Urban	0.005 (0.024)
Duty	-0.104 (0.021) ***
Interest	-0.035 (0.035)
Voted	-0.012 (0.065)
poly(Treatment, 2)1 × Independent	4.424 (2.077) *
poly(Treatment, 2)2 × Independent	-1.396 (2.085)
Num.Obs.	1761
R2	0.545
R2 Adj.	0.542
AIC	4628.5
BIC	4710.6
Log.Lik.	-2299.230
RMSE	0.89

予測値

コード

predictions(fit10,
            newdata = datagrid(Treatment   = 20:80,
                               Independent = 0:1)) %>%
  mutate(Independent = factor(Independent, levels = 0:1, 
                              labels = c("あり", "なし"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Independent), alpha = 0.25) +
  geom_line(aes(color = Independent), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "支持政党", fill = "支持政党") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit10, variables = "Treatment",
                newdata = datagrid(Independent = 0:1)) %>%
  mutate(Independent = factor(Independent, levels = 0:1, 
                              labels = c("あり", "なし"))) %>%
  ggplot(aes(x = Independent, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "支持政党", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit10, variables = "Independent",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "支持政党の有無の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.9 投票有無

コード

fit11 <- lm(Q16.1 ~ poly(Treatment, 2) * Voted + 
             Female + Age + Educ + Income + Urban + 
             Duty + Interest + Independent, 
   data = df)

modelsummary(fit11,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)

	(1)
(Intercept)	2.418 (0.159) ***
poly(Treatment, 2)1	45.016 (2.351) ***
poly(Treatment, 2)2	10.539 (2.486) ***
Voted	-0.024 (0.065)
Female	0.110 (0.045) *
Age	0.001 (0.001)
Educ	0.014 (0.025)
Income	-0.003 (0.004)
Urban	0.006 (0.024)
Duty	-0.105 (0.021) ***
Interest	-0.034 (0.035)
Independent	-0.021 (0.045)
poly(Treatment, 2)1 × Voted	-2.530 (2.615)
poly(Treatment, 2)2 × Voted	7.003 (2.729) *
Num.Obs.	1761
R2	0.546
R2 Adj.	0.543
AIC	4626.0
BIC	4708.1
Log.Lik.	-2297.996
RMSE	0.89

予測値

コード

predictions(fit11,
            newdata = datagrid(Treatment = 20:80,
                               Voted     = 0:1)) %>%
  mutate(Voted = factor(Voted, levels = 0:1, 
                        labels = c("棄権", "投票"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Voted), alpha = 0.25) +
  geom_line(aes(color = Voted), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率（%）", y = "投票率に対する評価 w/ 95% CI", 
       color = "投票有無", fill = "投票有無") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果（処置変数）

コード

marginaleffects(fit11, variables = "Treatment",
                newdata = datagrid(Voted = 0:1)) %>%
  mutate(Voted = factor(Voted, levels = 0:1, 
                        labels = c("棄権", "投票"))) %>%
  ggplot(aes(x = Voted, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "投票有無", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果（調整変数）

コード

marginaleffects(fit11, variables = "Voted",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数（= 投票率）の値（%）", 
       y = "投票有無の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

5 RDD

　40%が投票率評価の「高低」の分岐点?

コード

rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 40, data = df) %>%
  summary()

Call:
rdd::RDestimate(formula = Q16.1 ~ Treatment, data = df, cutpoint = 40)

Type:
sharp 

Estimates:
           Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|) 
LATE       3.192      265           0.4530    0.12171     3.722    0.0001977
Half-BW    1.596      118           0.3500    0.12496     2.801    0.0050972
Double-BW  6.384      487           0.3235    0.09812     3.297    0.0009785
              
LATE       ***
Half-BW    ** 
Double-BW  ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F-statistics:
           F      Num. DoF  Denom. DoF  p      
LATE       3.815  3         261         0.02115
Half-BW    2.536  2         115         0.16720
Double-BW  4.325  3         483         0.01010

コード

rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 40, data = df) %>%
  plot()

　ちなみにLATEが統計的に有意である閾値は40%と45%

コード

25:75 %>%
  enframe(name = "ID", value = "Cutpoint") %>%
  mutate(Fit = map(Cutpoint, ~rdd::RDestimate(Q16.1 ~ Treatment,
                                              cutpoint = .x,
                                              data = df)),
         Est = map_dbl(Fit, ~.x$est[1]),
         SE  = map_dbl(Fit, ~.x$se[1]),
         LL  = Est + qnorm(0.025) * SE,
         UL  = Est + qnorm(0.975) * SE,
         p   = map_dbl(Fit, ~.x$p[1]),
         Sig = if_else(p < 0.05, "Sig", "Insig")) %>%
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_pointrange(aes(x = Cutpoint, y = Est, ymin = LL, ymax = UL,
                      color = Sig)) +
  labs(y = "LATE w/ 95% CI") +
  scale_color_manual(values = c("Sig" = "black", "Insig" = "gray80")) +
  guides(color = "none") +
  theme_bw()

　むろん、この40%になにかの根拠があるわけではない。むしろ、change-in-slopeモデルで分析した場合、50%という結果も。

コード

pacman::p_load(cpop)

cpop(df2$Outcome, df2$Treatment)

cpop analysis with n = 61 and penalty (beta)  = 8.221748

1  changepoint detected at x = 
 50
fitted values : 
  x0       y0 x1       y1   gradient  intercept       RSS
1 20 1.113460 50 1.424349 0.01036297  0.9062005 0.4664509
2 50 1.424349 80 4.110136 0.08952625 -3.0519633 0.6733183

overall RSS = 1.139769
cost = 85.89115

コード

cpop(df2$Outcome, df2$Treatment) %>%
  plot()

　閾値を50%にしたRDDの場合、LATEは統計的に有意ではない。

コード

rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 50, data = df) %>%
  summary()

Call:
rdd::RDestimate(formula = Q16.1 ~ Treatment, data = df, cutpoint = 50)

Type:
sharp 

Estimates:
           Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|)   
LATE       3.36       234           -0.3896   0.3493      -1.115   0.2647     
Half-BW    1.68        82           -0.3333   0.2348      -1.420   0.1556     
Double-BW  6.72       459           -0.2948   0.2000      -1.474   0.1404     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F-statistics:
           F      Num. DoF  Denom. DoF  p      
LATE       2.118  3         230         0.19747
Half-BW    1.359  2          79         0.52593
Double-BW  4.270  3         455         0.01093

コード

rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 50, data = df) %>%
  plot()